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Abstract. This paper contains construction and analysis a finite element approx- 
imation for convection dominated diffusion problems with full coefficient matrix on 
general simplicial partitions in lR d , d > 2. This construction is quite close to the 
scheme of Xu and Zikatanov [28) where a diagonal coefficient matrix has been con- 
sidered. The scheme is of the class of exponentially fitted methods that does not use 
upwind or checking the flow direction. It is stable for sufficiently small discretiza- 
tion step-size assuming that the boundary value problem for the convection-diffusion 
equation is uniquely solvable. Further, it is shown that, under certain conditions 
on the mesh the scheme is monotone. Convergence of first order is derived under 
minimal smoothness of the solution. 



1. Introduction 

We consider the following convection-diffusion-reaction problem: Find u = u(x) 
such that 
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Here Q is a bounded polygonal domain in H d , d — 2,3, D — D(x) is d x d symmetric, 
bounded and uniformly positive definite matrix in Q, b* = (pi(x), . . . , bd(x)) is a 
given vector function, n is the unit outer vector normal to dQ, and / is a given source 
function. We have also used the notation Vw for the gradient of a scalar function u 
and V • b for the divergence of a vector function b in H d . The boundary of Q, dfl, is 
split into Dirichlet, T^, and Neumann, T^, parts. Further, the Neumann boundary is 
divided into two parts: T N = T™ U T^ 1 *, where T™ = {x G T N : n(x) ■ b(x) > 0} and 
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r^y* = {x G T N : n(x) • b(x) < 0}. We assume that T D has positive surface measure. 
The case D(x) = el, where / is the identity matrix in IR d and e > is a small 
parameter, corresponds to the important and difficult class of isotropic singularly 
perturbed convection-diffusion problems. 

Various generalizations have wide practical applications. For example, 7« could be 
replaced by nonlinear reaction term j(u) or the linear convective flux bu could be 
replaced by a nonlinear advection flux b(u). Finally, u could be a vector function 
describing the concentration of various chemicals or biological components so that 
(jl.ip is a system of equations coupled through the absorption/reaction term. Now 
7 is a matrix that models the chemical reactions or the biological interaction of 
the components. All these cases give rise to mathematical problems of convection 
dominated processes with possibly anisotropic diffusion. 

Our study of numerical method for solving (j 1.1ft is motivated by the fact that the 
above problem is the simplest model of transport and dispersion of a passive contam- 
inant in porous media. If the pressure p(x) in the aquifer is known (or already has 
been computed by solving a standard diffusion problem) then the pressure gradient 
forces the ground water to flow. The transport of a contaminant dissolved in the 
water, is described by the dispersion-reaction equation fll.ip . where u(x) represents 
the contaminant concentration, b = v = AVp is the Darcy velocity (up to a sign), 
A is the permeability of the porous media, 7 is the bio-degradation/absorption rate, 
and D(x) is the diffusion- dispersion matrix given by 

(1.2) D{x) = k d I + fctbb*/|t>| + fc,(|b|J - bb*/|b|). 

Here kd, k t , and ki are coefficients of diffusion, transverse dispersion, and longitudinal 
dispersions, respectively (cf. [in])- In dispersive underground flows k t > k[ which 
implies that D(x) is positive definite matrix, but possibly ill-conditioned. This prob- 
lem exhibits all difficulties associated with this class: monotone solutions that are 
highly localized due to internal and boundary layers, material heterogeneities and 
orthotropy, complex geometry, etc. 

Among the deficiencies of the standard finite element, finite volume, and finite 
difference approximations are loss of monotonicity, so that the numerical solution 
often exhibits non-physical oscillations, loss of solvability of the resulting algebraic 
problem, poor local resolution, fast dissipation of the energy, etc. A. A. Samarskii 
was one of the first to encounter the difficulties that arise in the numerical solution 
of such problems. In the early 60-ies A. A. Samarskii addressed most of the issues for 
one-dimensional problems that resulted in a new scheme described in his monograph 
[2UI Chapter 4]. 

In the past 40 years many special approximation techniques have been developed 
for multidimensional problems, for structured and unstructured grids, for general 
second order elliptic operators, etc. These techniques include monotone and upwind 
finite difference, finite volume, and finite element methods (e.g. [2], [8], [TTJ, [20], [2Tj . 
[22] . |23j . and (27]), streamline diffusion stabilization of the finite element method 
(e.g. [Tj, [3J, [5J, [13], and [33]), and special functional spaces setting (e.g. [7J and 
[24]). For more information regarding numerical methods and analytical techniques in 
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solving and studying convection-diffusion equations, especially convection dominated 
problems, we refer to the monograph of Ross, Stynes and Tobiska [19] . 

On a continuous level many convection-diffusion satisfy maximum principle. This is 
a desirable property of the solution of the resulting discrete problem as well. Scheme 
that satisfies maximum principle is often called monotone scheme. Among the sev- 
eral aforementioned schemes, upwind schemes are often monotone provided that the 
coefficient matrix D(x) is diagonal. In the recent works [2Tj, [22J Samarskii and his 
co-workers were able to derive monotone schemes on rectangular meshes for the prob- 
lem (II. ip when D(x) is a full matrix. These schemes are second order accurate on 
uniform meshes and solution in C 3 . 

The idea of construction of monotone schemes for singularly perturbed convection- 
diffusion problem goes back to the work by Scharfetter and Gummel [25], where the 
monotonicity has been a very desired property in numerical semiconductor device 
modeling. Exponentially fitted scheme for a general convection-diffusion problem 
with diagonal matrix D(x) on an arbitrary simplicial mesh was derived and studied 
in [12], [I], [5], [IS]) [IB], [17] . As we mentioned, the technique we use here is a 
generalization of [28] for general meshes and is also related to the methods proposed 
in [30] (see also [3]] for simple diode simulation using such scheme). Under mild 
conditions on the mesh (the partition has to satisfy certain angle condition) it has 
been shown that the scheme is monotone. Further, in [28] it was proved that the 
scheme converges with first order provided that the solution u G W 1,p and the flux 
D(x)V« + b(x)u G (W l,p ) d for p > d. Note, that these are very mild conditions on 
smoothness of the solution of problem ( II. ip . 

The aim of this note is to construct an exponentially fitted finite element ap- 
proximation of (II .ip on general simplicial partitions, for symmetric positive definite 
matrices D(x), and for arbitrary vector- functions b. The proposed scheme is a gener- 
alization of the discretization derived in [28] , [30] for problems with diagonal matrices 
D(x). Important role in the construction and the analysis plays the expansion of a 
constant over each element vector-flux using the lowest order Nedelec basis for sim- 
plicial finite elements. This allows to present the bilinear form through differences of 
the vertex values of the test functions and the exponentially weighted local solution. 
This representation ensures the consistence of the method. 

The scheme has several interesting features. It is a finite element scheme with a 
standard variational formulation (but with a modified bilinear form); it does not use 
explicitly the standard upwind techniques, such as checking the flow directions; it can 
be applied to very general unstructured grid in any spatial dimension. It would be 
difficult to expect that in such generality the scheme will be monotone. Nevertheless, 
we were able to find conditions that involve the geometry of the finite elements in a 
metric associated with D so that the scheme is monotone. Further, for sufficiently 
small step-size of the finite element partition we prove existence and uniqueness of 
the solution of the discrete problem by using a fundamental result of Schatz [25] . 

The paper is organized as follows. In Section |2J we introduce the necessary notation 
for Sobolev spaces, finite element partition and the discrete space. Section [3] contains 
the main results of the paper. In Subsection 13.11 we present the rationale used in 
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derivation of exponentially fitting finite element scheme. An important concept here 
is the edge based interpolation of the total flux that uses an ordinary differential equa- 
tion along the edge. In Subsection 13.21 we present the scheme itself as a consequence 
of this special interpolation. The main result here is contained in Lemma 13.11 where 
certain properties of the discrete bilinear form are obtained. Finally, in Subsection 
13.31 we prove the stability of the scheme for sufficiently small mesh-size and derive an 
estimate for the error under minimal smoothness of the solution. 



2. Notation 

In this section, we introduce the necessary notation and describe some basic prop- 
erties of finite element partitions and finite element spaces. 

We denote by L P (K), 1 < p < oo the space of p-integrable real- valued functions 
over K C Q (with the usual modification for p = oo), by (•, -)k and \\-\\k, respectively, 
the inner product and the norm in L 2 (K). Further | • \i )P) k an d || • respectively 
denote the semi-norm and norm of the Sobolev space W 1,P (K). For p = 2 we use 
H X {K) := W 1,2 (K) and if K = Q often we suppress the index K so that (•, -)q := (•, •) 
and || • ||ft := || • ||, and || • \\i t n := || • ||i. Further, we use the Hilbert space 

H 1 D (Q) = {veH 1 (Q): v\ Vo = 0}. 

We introduce the bilinear form a(-, •) defined on Hp(Q) x iJ^(fi): 



(2.1) a(w, u) := (DVm + bw, Vv) + (7M, - / b ■ n uvds. 

J 1 N 

Then (II -ip has the following weak form: Find u G Hjj(Q) such that 

(2.2) a{u,v) = F(v) := (f,v) + f gvds for all v e Hp (Q) . 

Jr™ 

Further in the paper we assume that the following inf-sup condition is valid: there is 
a constant cq > 0, such that 

(2.3) sup > co\\w\\ h Vw6 4(0). 

veH^(n) IMIi 

We shall also assume that the bilinear form a(w,v) is bounded on H}j(Q) and the 
linear form F(v) is continuous in Hj~,(Q). Then the above problem has unique solution 
(cf. [12]). 

Remark 2.1. A sufficient condition for (12. 3p and continuity of a(u, v ) and F(v) are, 
for example, 7(0;) + 0.5V • b(a;) > for all x G Q, boundedness of the coefficients 
D(x), b(x), and 7(2) in Q. 

Let Th be a family of simplicial finite element triangulations of Q that are shape 
regular and satisfy the usual conditions (see [91 Chapter 2]). For simplicity of the 
exposition, we assume that the triangulation covers Q exactly. Associated with each 
Th, let Vh C Hp(Q) be the finite element space of piece- wise linear functions. By 
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v i G Vh we denote the standard finite element Lagrange interpolant which assumes 
the values of v G C° at the vertexes in the partition Th- 

Given T G Th, we introduce the following notation. By qj, j — 1, . . . , 4 we denote 
the vertices of T, E is the edge connecting two vertices qi and qj, 5e4> = 4>(.<li) ~ ^(Qj) 
for any continuous function <fi on E, and t e = 5 E x = q,i — qj is a directional vector of 
E (not assumed unitary). 

3. Exponential fitting scheme for general convection-diffusion 

equations 

3.1. Preliminaries. Introduce a notation for the scaled flux 

(3.1) J(u) =Vu + (3(x)u, /3 = D- 1 b. 

We will assume further that J{u) G [W l > p {tt)] d , p > n, D, D" 1 G [W 1 ' 00 (tt)] dxd and 
b G W 1,oc (Q). These assumptions on the coefficients smoothness can be relaxed to 
hold element-wise (i.e. for each T G Th) and the considerations below will still hold 
with changes of some of the norms used in the error estimate to be taken element- wise 
as well. 

The basic idea which we use in the construction of the exponentially fitted scheme 
is to approximate the flux vector J(u) with a constant vector field Jt{u) on each 
element T of the partition Th- Apparently, if Jt{u) is a constant on each simplex, 
then we can expand it using the Nedelec basis as follows: 

J T = ^Jt- t e <Pe(x)- 

Here fE are the Nedelec basis functions, which in terms of the barycentric coordinates 
Xi are given by 

y? E := XiVXj - XjVXi, E = (q h qj). 
The goal then is to write out Jt(u) ■ r E in terms of u, for all edges E and thus 
determine the approximation. To find the moments of the tangential flux, we use the 
same technique as in |2H]. Let u G H^(Q) nC°(0). Consider an edge e C dT. Taking 
the Euclidean inner product with te we obtain 

(Vu ■ r E ) + (f3 ■ t e )u = (J(u) ■ t e ). 

A change of variables in this ordinary differential equation then gives: 

(3.2) e -^dE(e^u) = T ^(J(u)-T E ), where «9 E ^ E = tK(P • ^e) 

I t e| |t e | 

and 8ev := Vi>-te/|te| is the directional derivative along the edge E. After integration 
over E we obtain that 

S E (e^u) = r^-r [ e^{J{u) ■ r E )ds. 
Fe| Je 

Let He(/3) be the harmonic average of e^ E over E defined as follows: 



(3.3) U E ((3) 



\ t e\ Je 



e^ds 
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The constant approximation J? is then obtained by using the mean value theorem 
J* ■ r E J E e^E ds — J E J- r E e^E ds, and the definition then is 

J T (u) ■ r E := J* • r E = H E (f3)5 E (e^u). 

3.2. Discrete problem. We now have all the ingredients needed to define the dis- 
crete approximation to (12. 2p . Based on the above considerations, we shall define two 
approximate bilinear forms. The first one is used in the formulation of the discrete 
problem and the second is used in an intermediate step needed to prove the error 
estimate. 

On a fixed element T C 7/,, we first introduce 

(3.4) a hjT (u h ,v h ) = ^(D)n E ((3)5 E (e^u h )5 E v h , 

EcdT 

where 

uil(D) = - J DV\ ■ V\j dx, E = (q u qj ). 

Note, that u^(D) give the element stiffness matrix for the diffusion part of the dif- 
ferential equation, —V • (D(x)Vu). 

Next, we use the expansion via the Nedelec basis functions, to define 

(3.5) b h>T {u h ,v h ) = U E {P)8 E {e^u h ) f D Ve ■ Vv h dx. 

EcdT ^ T 

The global bilinear form is then obtained by summing over all elements of the trian- 
gulation the local forms (13. 4p and adding the contributions from the boundary T™ 1 , 
as follows: 

(3.6) a h (u h ,v h ) = ^2 a h , T (u h ,v h ) + / ^u h v h dx - / b-nu h v h ds. 

TcT h Jn ecv^ Je 

Finally, the finite element approximation of the problem (jl.ip reads as follows: Find 
Uh G Vh such that 

(3.7) a h (u h ,v h ) = F(v h ), for all v h E V h . 

The following lemma is the main tool used in the analysis of the above scheme. 

Lemma 3.1. The following relations hold for any Vh G V^: 

1. If we C(T) then 

(3.8) b h , T (wj,v h )= V VMS. f e ^J( w ). TE ds\ [ Dip E .Vv h ; 

EcdT L \ t e\ Je J Jt 

2. If Jrp is a constant vector on T, then for any Vh G 14 

(3.9) ^2 Jt-t e Dip E ■ Vv h dx = ^ u e( d ) j t • T E 5 E v h ; 

EcdT ^ T EcdT 
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3. Ifw G C(T) and J{w) G [H /1,p (T)] n , p > d, then the following inequality holds 
for every Vh G T4 and T G Th 



(3.10) 



\a T (w,v h ) - a h ,T(wi,v h )\ < Ch\J(w)\ liPiT \\v h \\ liT . 



where 



a T (w,v h 



(DVw + hw) ■ Vvh dx. 



T 



Proof. The proof of 1. follows directly from the derivation. 

The proof of 2. can be done as follows: Consider $ := Jt ■ x for x G T (here x is 
the vector of coordinates in ]R d ). It is obvious that $ is linear and that V$ = Jt and 
hence V$-te = Jt-te- We now use the fact that the Nedelec canonical interpolation 
Jt = J^ect^t ■ t e)^Pe satisfies the commutativity property ITa/V^ = V$/ = V$. 
Therefore, 



^2 Jt-t e Dip E ■ Vv h dx = 

JT JT 



D 



EcdT 



2J Jt ■ te¥>e 

EcdT 



Vvh dx 



D[U M J T ]-Vv h dx 



[ D[U M J T ] ■ Vv h dx= £>[n^v$] ■ Vv h dx 
Jt Jt 

/ 7J[W] ■ Vv h dx = Ue(D)J t ■ te5ev h , 

^ T EcdT 



and the proof of 2. is complete 

To prove 3. we use (13.81) and split the difference in the following way 



(3.11) 
where 

and 



a T (w,v h ) - a h)T {wi,v h ) = £i{J{w),v h ) + £ 2 (J(w),v h ) 
£i(J(w),v h ) = a T (w,v h ) - b hiT {wi,v h ), 
£ 2 {J(w),v h ) = b h)T {wi,v h ) - a hjT (wi,v h ). 



From the relations and item 1. we can expand the forms a T (w,Vh), CLh^iwi^h) and 
bh,T(uii,Vh) to get that 



£t{J{w),v h ) 



(3.12) 



J{w) • Vvhdx 



E 

EcdT 



u E {P) 



f e^J(w)-T E ds 
Je 



D<p E ■ Vvh dx 
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£ 2 (J(w),v h ) 



Ec&r 



|T E 



/ e^J( w ) ■ r E ds 
Je 



Dip E ■ Vvh dx 



EcdT 



T E 



[ e^J(w)-r E ds 
Je 



A change of variable from the standard reference element T (of unit size) to T: 
x = Bx + 6q and w(x) 



a T (w,v h 



ah,T{w,v h ) 



bT,h{w,v h ) 



w(x) gives the following scaled bilinear forms 
= |det£| J (B~ 1 dJ(w) ■ Vv h ) dx 



Ecf 

\detB\ E 



Ue(P) 



EcT 



X 



FeI 



Z^{J{w) -T E )ds 



SeVh 



J^B-'Difv ■ Wv h ) dx. 



By our assumptions on the smoothness of J(w), the corresponding error functionals 
£i(J(w),Vh), i — 1,2, can be appropriately bounded: 



(3.14) 



£i(J(w),v h ) < Ci\\J(w)\\ t \\v h \\ h f 



T- 



where Ci might depend on D, but do not depend on (3. By the Sobolev inequality 
we have that 

||JHIIo,oo,T <C|I^HIIl >P> f, P >d. 

We observe that from (13.81) and (13.91) . it follows that £i(J(w),Vh) = if J{w) is a 
constant on T. By applying the Bramble-Hilbert Lemma on T, and scaling back to 
T we obtain the desired result: 



(3.15) 



\£i(J(w),v h )\ < Ch\J(w)\ 1)P)T \v h \ 



1,T> 



1,2. 



□ 



3.3. Solvability of the discrete problem and error estimate. In this paragraph 
we state two lemmas related to the solvability of the problem and then a result related 
to the error bound. The first result, the proof of which follows straightforward from 
the definition, is related to the monotonicity of the scheme (i.e. discrete maximum 
principle). This amounts to a condition on the geometry of the mesh associated with 
the matrix D. 

Lemma 3.2. The stiffness matrix corresponding to the bilinear form ( [<?. 7| ) is an Ad- 
matrix for any continuous function f3 if and only if the following inequality holds for 
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all edges E in the triangulation 

(3.16) $>^° 

One may check out easily that if D is a constant matrix and d = 2 (two spatial 
dimensions) this is equivalent to the statement that the triangulation is Delaunay in 
the metric introduced by D. Namely, instead of Euclidean inner product b • n of the 
vectors b and n in H d we need to use the inner product Dh ■ n (recall that D is a 
symmetric and positive definite matrix). In this case the global stiffness matrix of 
the finite element system is nonsingular and therefore the scheme (13. 7p has unique 
solution. 

Next result is about solvability of the discrete problem for sufficiently small char- 
acteristic mesh size h. Let us first consider an auxiliary discrete problem with a(-, ■) 
in place of ah(-, •) in (13. 7p . The latter problem is solvable, and a convincing (but 
not rigorous) argument to prove this claim is that the convection term is one order 
lower than the diffusion term and hence, decreasing h will make the diffusion term 
dominating and the problem weakly coercive. Some more detailed considerations and 
a rigorous arguments can be found in Schatz [26J or Xu 



Lemma 3.3. For sufficiently small h the following inf-sup condition holds 
/o i7\ a h (w h ,v h ) 

(3.17) sup — | — |i > cillwfeU^n vw h eV h 

v h ev h \\Vh\\i,a 

with a constant c\ > independent of mesh-size h. 

Proof. As we have pointed out, when the original bilinear form a(-, •) is used in (13. 7p . 
the discrete problem is uniquely solvable (for sufficiently small h). Hence, there exists 
a constant C2 such that 

(3.18) sup -j: — I, >c 2 \\w h \\i i si, Vw h eV h . 

v h ev h \\Vh\\i,n 

Let Vh,Wh € Vh- Then obviously 

a h (wh, v h ) = a(w h , v h ) + [a h (w h , v h ) - a(w h , v h )\ . 

The first term is estimated using the condition (I3.18p . To estimate the second term 
we use (I3.10p from lemma I3.1[ sum up over all T and apply the Schwarz inequality 
to obtain that 

\a(w h ,v h ) - a h {w h ,v h )\ < Ch I ^ |J r (wh)\i, p ,T \ \\ v h\\i,n- 
Observing that \wh\2,T = for any uih G Vh, T G Th we get 

Summing over all the elements of the partition we have 

(3.19) \a(w h ,v h ) - a h (w h ,v h )\ < Ch max ||j9||i j0 o/r||wA||i,nlKlli,n, 

TeTh 



max||/3||i i00 , T 



2 P,T 
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and for h satisfying 

h < ho = C 

the discrete problem has a unique solution. □ 
As a consequence of Lemma I37TI Lemma [3731 we get the following convergence result. 

Theorem 3.4. Let u be the solution of the problem ( QQ| ). Assume that for all T G 7h, 

D G (W 1 ' 00 ^))^^ /3 G [W rl '°°(T)] d ; it G W^T), 7 G C(T) ; and J(w) ee Vu + 
(3(x)u G (H /1,p (T)) d , p > d. Then for sufficiently small h, the following estimate 
holds: 

(3.20) Utij - UfcH^n < C7i i J2 \ J ( u )\i,p,t + Mi 

Remark 3.5. There are also other possibilities for expressing the flux J(w). For ex- 
ample, instead of the flux ( 13.1 j) one can write 

D(x)Vu + bu = D(a;)a(x)~ 1 (a(a;)V(u) + a(x)D" 1 (a;)bu)) 
:= D(ar) (a(x)V(u) + /3m)) , 

where 

/3 = a(x)D- 1 (x)b, = D(x)a(:c) _1 

with a(x) a suitable positive function (or a positive diagonal matrix). Then define 

J(u) = a(x)V(u) + flu. 

For such a choice of J(u) the derivation and the analysis of an exponentially fitting 
scheme are essentially the same with some changes occurring in the harmonic averages 
used to define the discrete problem. 
For example, one may choose 

a(x) = (X min (D(x)) + X max (D(x))/2, 

where X m i n (D(x)) and X max (D(x)) are the minimum and maximum eigenvalues of 
D(x). For such choice D^ 1 = aD^ 1 is better conditioned. For example, problem 
(HZX), (HD with data such as k d = 0.0001, k t = 21, and h = 2.1, used in jg], might 
require such modification. 

Remark 3.6. As we have pointed out in the introduction, in many cases D(x) takes 
the form (jl.2p . Then introducing the orthogonal projection 7Tb = bb'/|b| 2 along the 
vector b(x) we can rewrite D(x) in the form 

D(x) = k d I + k t \b\ii h + ki\b\(I - n h ). 

Now one easily finds that D^h = (kd + /^|b|) -1 b, i.e. the evaluation of D~ x h is just 
a multiplication of b by a scalar. 
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